function[X,Y,Z,X1,Y1] = rob_regress(X,Y,Z,X1,Y1),
     [rs,rstat] = robustfit(X,Y);
     [rs1,rstat1] = robustfit(X,Z);
     [rs2,rstat2] = robustfit(X,X1);
     residuals = abs(rstat.resid);
     residuals1 = abs(rstat1.resid);
     residuals2 = abs(rstat2.resid);
     LLL = residuals == max(residuals) | residuals1 == max(residuals1); %Just remove the most noisy point
     X(LLL) = [];
     Y(LLL) = [];
     Z(LLL) = [];
     X1(LLL)= [];
     Y1(LLL) = [];
    % Robust outlier rejection  --- END